We consider the Hamiltonian
$$ H(r, \theta, p_r, p_\theta) = p_r\, \mu(r) + \Vert p \Vert_{g} $$
where
$$ \mu(r) = \delta \sin 2r, \quad \Vert p \Vert_{g} = \sqrt{ p_r^2 + \frac{p_\theta^2}{m_\lambda^2(r)} }, \quad m_\lambda^2(r) = \frac{\sin^2 r}{1 - \lambda \sin^2 r} $$
and where $\delta$ and $\lambda$ are constants. We will fix $\delta > 1$ (strong Zermelo case) and $\lambda = 4/5$. Along the geodesics, we have $H+p^0 = 0$. The parameter $p^0$ is constant equal to $-1$ (hyperbolic), $0$ (abnormal) or $1$ (elliptic).
Remark. We can parameterize the geodesics by the norm of the initial convector, setting $\Vert{p_0}\Vert_g = 1$. This amounts to parameterize by the initial angle $\alpha_0$: $$ p_r = \sin \alpha_0, \quad p_\theta = m_\lambda(r) \cos \alpha_0. $$ In that case, the hyperbolic geodeics are given by $$ p_r\, \mu(r) + 1 = \sin(\alpha_0)\, \delta \sin(2r) + 1 > 0. $$
import sys
sys.path.insert(1, '../')
#
import math
import numpy as np # scientific computing tools
import wrappers # for compilation of Fortran Hamiltonian codes
import geometry2d.database
import geometry2d.conjugate
import geometry2d.errors
import geometry2d.geodesic
import geometry2d.plottings
import geometry2d.problem
import geometry2d.splitting
import geometry2d.utils
import geometry2d.wavefront
import matplotlib.pyplot as plt # for plots
%matplotlib inline
#
color_hyperbolic = 'red'
color_elliptic = 'blue'
color_abnormal = 'green'
color_cut_locus = 'black'
color_conjugate = 'magenta'
#
color_strong_domain_2d = 'lightblue'
color_strong_domain_3d = 'blue'
# print Fortran code of the Hamiltonian: print hfun.f90 file
with open('hfun.f90', 'r') as f:
print(f.read())
! Lindblad
subroutine hfun(x, p, d, h)
double precision, intent(in) :: x(2), p(2), d
double precision, intent(out) :: h
! local variables
double precision :: r, th, pr, pth
double precision :: l, m2, mu
r = x(1)
th = x(2)
pr = p(1)
pth = p(2)
l = 4d0/5d0
m2 = sin(r)**2/(1d0-l*sin(r)**2)
mu = d*sin(2*r)
h = mu*pr + sqrt(pr**2 + pth**2/m2)
end subroutine hfun
def initialize(d, r0):
# Parameters
#d = -1.25 # strong Zermelo case
#r0 = np.pi/2.0 #+ 3.0*np.pi/8.0 # initial azimuth
θ0 = 0.0 # initial latitude
t0 = 0.0 # initial time
name = 'lindblad_strong' # name of the problem
# Initialize data
data_file = 'data_lindblad_strong.json'
restart = False # restart or not the computations
data = geometry2d.database.Data({'name': name,
't0': t0,
'r0': r0,
'θ0': θ0,
'd': d}, data_file, restart)
# Initial point
q0 = np.array([r0, θ0])
# Hamiltonian and derivatives up to order 3
H = wrappers.hamiltonian(d, compile=False, display=False)
# The Riemannian metric associated to the Zermelo problem
# g = g1 dr^2 + g2 dθ^2
def g(q):
λ = 4.0/5.0
r = q[0]
g1 = 1.0
g2 = np.sin(r)**2/(1.0-λ*np.sin(r)**2)
return g1, g2
# problem
prob = geometry2d.problem.GeometryProblem2D(name, H, g, t0, q0, data, steps_for_geodesics=200)
return prob
# function to plot the domain of strong current between two angles r1 and r2
def plot_2d_domain(fig, r1, r2):
# from r to phi
def phi(r):
return r - np.pi/2.0
ph1 = phi(r1)
ph2 = phi(r2)
# plot the domain: a 2d surface
ax = fig.axes[0]
x = np.linspace(-np.pi, np.pi, 100)
y1 = np.ones_like(x)*ph1
y2 = np.ones_like(x)*ph2
ax.fill_between(x, y1, y2, color=color_strong_domain_2d, alpha=0.5)
return fig
def plot_3d_domain(fig, r1, r2, *, elevation, azimuth):
#
color = color_strong_domain_3d
alpha = 0.5
# from r to phi
def phi(r):
return r - np.pi/2.0
φ1 = phi(r1)
φ2 = phi(r2)
# plot the domain: a 2d surface between the parallel of latitude φ1 and φ2 on the sphere
ax = fig.axes[0]
#
N = 100
u = np.linspace(0, 2 * np.pi, N)
v = np.linspace(φ1, φ2, N)
x = np.outer(np.cos(u), np.cos(v))
y = np.outer(np.sin(u), np.cos(v))
z = np.outer(np.ones(np.size(u)), np.sin(v))
#
x_front = np.zeros((N, N))
y_front = np.zeros((N, N))
z_front = np.zeros((N, N))
x_back = np.zeros((N, N))
y_back = np.zeros((N, N))
z_back = np.zeros((N, N))
cam = geometry2d.plottings.get_cam(elevation, azimuth, geometry2d.plottings.dist__)
gap = 1e-1
for i in range(N):
for j in range(N):
u = np.array([x[i,j], y[i,j], z[i,j]])
v = np.array([cam[0], cam[1], cam[2]])
ps = np.dot(u,v)/(np.linalg.norm(u)*np.linalg.norm(v))
#ps = x[i,j]*cam[0]+y[i,j]*cam[1]+z[i,j]*cam[2]
if ps >= -gap:
x_back[i,j] = x[i,j]
y_back[i,j] = y[i,j]
z_back[i,j] = z[i,j]
x_front[i,j] = np.nan
y_front[i,j] = np.nan
z_front[i,j] = np.nan
if ps <= gap:
x_front[i,j] = x[i,j]
y_front[i,j] = y[i,j]
z_front[i,j] = z[i,j]
x_back[i,j] = np.nan
y_back[i,j] = np.nan
z_back[i,j] = np.nan
#
ax.plot_surface(x_front, y_front, z_front, color=color, alpha=1) #, zorder=z_order_surface)
ax.plot_surface(x_back, y_back, z_back, color=color, alpha=alpha) #, zorder=-z_order_surface)
return fig
# we parameterize ||p(0)|| = 1 by the angle α0
# the norm being given by the metric g: p(0) = (sin(α0)*g1, cos(α0)*g2) = (sin(α0), cos(α0)*g2), since g1 = 1
# get the values of α0 for which h = pr mu(r) + 1 > 0
#
# first get the values of α for which h = pr mu(r) + 1 = 0
def mu(r, d):
return d*np.sin(2*r)
def pr(α):
return np.sin(α)
def h(r, α, d):
return pr(α)*mu(r, d) + 1.0
def α0(r, d):
return np.arcsin(-1.0/mu(r, d))
# get the values of α0 for which h = pr mu(r) + 1 > 0
def α_hyperbolique_1(r, d):
m = mu(r, d)
if abs(m) > 1.0:
α = α0(r, d) # between -pi/2 and pi/2
if α > 0.0:
return np.pi - α
else:
return α
else:
return 0.0
def α_hyperbolique_2(r, d):
m = mu(r, d)
if abs(m) > 1.0:
α = α0(r, d) # between -pi/2 and pi/2
if α > 0.0:
return 2.0*np.pi + α
else:
return np.pi - α
else:
return 2.0*np.pi
def α_hyperbolique_span(r, d, *, N=20):
α1 = α_hyperbolique_1(r, d)
print('α1 = ', α1)
α2 = α_hyperbolique_2(r, d)
print('α2 = ', α2)
return np.linspace(α1, α2, N)
# get the values of α0 for which h = pr mu(r) + 1 < 0
def α_elliptique_1(r, d):
m = mu(r, d)
if abs(m) > 1.0:
α = α0(r, d) # between -pi/2 and pi/2
if α > 0.0:
return α
else:
return np.pi - α
else:
# raise an error: there is no elliptic domain
raise geometry2d.errors.ArgumentValueError('There is no elliptic domain for this value of d')
def α_elliptique_2(r, d):
m = mu(r, d)
if abs(m) > 1.0:
α = α0(r, d) # between -pi/2 and pi/2
if α > 0.0:
return np.pi - α
else:
return α + 2.0*np.pi
else:
# raise an error: there is no elliptic domain
raise geometry2d.errors.ArgumentValueError('There is no elliptic domain for this value of d')
def α_elliptique_span(r, d, *, N=20):
α1 = α_elliptique_1(r, d)
print('α1 = ', α1)
α2 = α_elliptique_2(r, d)
print('α2 = ', α2)
return np.linspace(α1, α2, N)
def make_plots(r0, d, prob, *, cameras=[{'azimuth': 160, 'elevation': -10}]):
# North hemisphere
rp1 = math.asin(1/np.abs(d))/2.0
rp2 = (math.pi - math.asin(1/np.abs(d))) / 2.0
# South hemisphere
rm1 = (-math.asin(1/np.abs(d)) % (2*np.pi) )/2.0
rm2 = (math.pi + math.asin(1/np.abs(d))) / 2.0
# Geodesics
geodesic = geometry2d.geodesic.Geodesic(prob)
#
αspan_hyperbolique = α_hyperbolique_span(r0, d, N=31)
#
fig2d = geodesic.plot(alphas=αspan_hyperbolique, tf=4, length=1.0,
view=geometry2d.plottings.Coords.PLANE, figsize=(4,4), color=color_hyperbolic)
fig2d = plot_2d_domain(fig2d, rp1, rp2)
fig2d = plot_2d_domain(fig2d, rm1, rm2)
# cameras is a list of dictionaries containing azimuth and elevetion angles
figs3d = []
for cam in cameras:
azimuth = cam['azimuth']
elevation = cam['elevation']
fig3d = geodesic.plot(alphas=αspan_hyperbolique, tf=4, length=1.0,
view=geometry2d.plottings.Coords.SPHERE, azimuth=azimuth, elevation=elevation, figsize=(3,3),
color=color_hyperbolic)
fig3d = plot_3d_domain(fig3d, rp1, rp2, azimuth=azimuth, elevation=elevation)
fig3d = plot_3d_domain(fig3d, rm1, rm2, azimuth=azimuth, elevation=elevation)
figs3d.append(fig3d)
return fig2d, figs3d
d = 1.25
r0 = np.pi/2.0
prob = initialize(d, r0)
fig2d, figs3d = make_plots(r0, d, prob)
α1 = 0.0 α2 = 6.283185307179586
<Figure size 640x480 with 0 Axes>
fig2d
figs3d[0]
d = 1.25
r0 = np.pi/2.0 + 2.0*np.pi/8.0
prob = initialize(d, r0)
cameras = [{'azimuth': 160, 'elevation': -10}, {'azimuth': 160, 'elevation': -10}]
fig2d, figs3d = make_plots(r0, d, prob, cameras=cameras)
# add Elliptic geodesics
αspan_elliptic = α_elliptique_span(r0, d, N=20)
geodesic = geometry2d.geodesic.Geodesic(prob)
fig3d_with_elliptic = geodesic.plot(alphas=αspan_elliptic, tf=np.pi, \
view=geometry2d.plottings.Coords.SPHERE, figure=figs3d[1], color=color_elliptic)
α1 = 2.214297435588181 α2 = 7.2104805251811985 α1 = 0.9272952180016123 α2 = 2.214297435588181
<Figure size 640x480 with 0 Axes>
fig2d
figs3d[0]
fig3d_with_elliptic
d = 1.25
r0 = np.pi/2.0 + 3.0*np.pi/8.0
prob = initialize(d, r0)
fig2d, figs3d = make_plots(r0, d, prob)
α1 = 0.0 α2 = 6.283185307179586
<Figure size 640x480 with 0 Axes>
fig2d
figs3d[0]
d = -1.25
r0 = np.pi/2.0
prob = initialize(d, r0)
cameras = [{'azimuth': 160, 'elevation': -10}, {'azimuth': 20, 'elevation': -25}]
fig2d, figs3d = make_plots(r0, d, prob, cameras=cameras)
α1 = 0.0 α2 = 6.283185307179586
<Figure size 640x480 with 0 Axes>
fig2d
figs3d[0]
figs3d[1]
d = -1.25
r0 = np.pi/2.0 + 2.0*np.pi/8.0
prob = initialize(d, r0)
cameras = [{'azimuth': 160, 'elevation': -10}, {'azimuth': 160, 'elevation': -40}]
fig2d, figs3d = make_plots(r0, d, prob, cameras=cameras)
α1 = -0.9272952180016123 α2 = 4.068887871591405
<Figure size 640x480 with 0 Axes>
fig2d
figs3d[0]
figs3d[1]
d = -1.25
r0 = np.pi/2.0 + 3.0*np.pi/8.0
prob = initialize(d, r0)
cameras = [{'azimuth': 160, 'elevation': -10}, {'azimuth': 160, 'elevation': -40}]
fig2d, figs3d = make_plots(r0, d, prob, cameras=cameras)
α1 = 0.0 α2 = 6.283185307179586
<Figure size 640x480 with 0 Axes>
fig2d
figs3d[0]
figs3d[1]